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In the framework of the paradigmatic prisoner's dilemma game, we investigate the evolutionary 
dynamics of social dilemmas in the presence of "cooperation facilitators" . In our model, cooperators 
and defectors interact as in the classic prisoner's dilemma, where selection favors defection. However, 
here the presence of a small number of cooperation facilitators enhances the fitness (reproductive 
potential) of cooperators, while it does not alter that of defectors. In a finite population of size 
N, the dynamics of the prisoner's dilemma with facilitators is characterized by the probability that 
cooperation takes over (fixation probability) and by the mean times to reach the absorbing states. 
These quantities are computed exactly and using Fokker-Planck equations. Our findings, corrob- 
orated by stochastic simulations, demonstrate that the influence of facilitators crucially depends 
on the difference between their density z and the game's cost-to-benefit ratio r. When z > r, the 
fixation of cooperators is likely in a large population and, under weak selection pressure, invasion 
and replacement of defection by cooperation is favored by selection if b(z — r)(l — z) > N~ , where 
< b < 1 is the cooperation payoff benefit. When z < r, the fixation probability of cooperators is 
exponentially enhanced by the presence of facilitators but defection is the dominating strategy. 
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I. INTRODUCTION 



Understanding the origin of cooperative behavior is a 
central issue in the life and behavioral sciences, and has 
recently been listed among the major scientific puzzles to 
be elucidated Jl). Evolutionary game theory (EGT) pro- 
vides the ideal framework to study the competition be- 
tween species and there is a long tradition of modeling the 
evolution of cooperation using evolutionary games ^ [| . 
In recent years, these processes have increasingly been 
investigated using the methods of statistical physics, see 
e.g. Q and references therein. In EGT, successful species 
spread at the expense of the others, and each individ- 
ual's reproductive potential (fitness) varies with the pop- 
ulation's composition that continuously changes in time. 
The interaction between the species is thus accounted for 
by a fitness-dependent (or "frequency-dependent" ) selec- 
tion pressure 0] , as observed in various experiments [Q . 
Quite intriguingly, in such a setting the optimization of 
the fitness at an individual level can result in the re- 
duction of the population overall fitness ^ ||. An in- 
fluential example of such a paradoxical behavior, is pro- 
vided by the celebrated prisoner's dilemma (PD) game 
that serves as a metaphor for social dilemmas. In fact, in 
the classic PD individual interest leads to defection, even 
though mutual cooperation would be socially more ben- 
eficial |^|, ||. While the PD is the paradigmatic model 
for the evolution of cooperation, its main prediction is 
at odds with the cooperative behavior that is commonly 
observed in experimental realizations [Q, ^). This has 
motivated an upsurge of research aiming to identify the 
possible mechanisms capable of promoting cooperation 
in biological and social systems Jq|. Notably, it has been 



proposed that cooperation can be promoted by kin and 
group selection [Q, as well as by conditional behavioral 
rules leading to direct or indirect reciprocity @, 0], It 
has also been found that local interactions may pro- 
mote cooperation in some social dilemmas . Further- 
more, it has been shown that cooperation is supported in 
games with voluntaryparticipation, or with punishment 
for non-cooperation |pl| . 

In this work, we investigate an alternative scenario for 
the spread of cooperation in social dilemmas: we con- 
sider the evolution of the prisoner's dilemma in a finite 
population comprising a small number of "cooperation 
facilitators" . The facilitators participate in the dynam- 
ics only by enhancing the reproductive potential of coop- 
erators, while they do not affect the fitness of defectors 
(see Sec. II below). To study the influence of cooperation 
facilitators on the prisoner's dilemma dynamics, the evo- 
lution is modeled in terms of a birth-death process and 
the fixation properties are studied analytically. In fact, 
it is well established that the evolutionary dynamics in 
finite populations is efficiently characterized by the prob- 
abilities of reaching the absorbing states, where the ex- 
tinction of one or more species and the fixation of another 
occur H 12 15 . Here, we are particularly interested in 
the probability that, from a given initial composition, the 
population eventually comprises only cooperators and a 
small fraction of facilitators, but no defectors ("cooper- 
ation fixation probability"). The mean times for these 
events (mean fixation times) are also studied and our 
results are checked against stochastic simulations. This 
approach allows us to (i) discuss how demographic fluc- 
tuations alter the mean field predictions of the classic 
replicator equations 0, and (ii) thoroughly analyze the 
circumstances under which facilitators and selection favor 
a single cooperator invading and replacing a population 
of defectors. 
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This paper is organized as follows: The PD with co- 
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operation facilitators is introduced in the next section, 
where some of its properties are discussed. In Section 
III the dynamics with the Fermi process is characterized 
by the fixation probability (Sec. III. A) and the mean fix- 
ation times (Sec. III.B). The dynamics with the Moran 
process is studied in Section IV, while we summarize our 
findings and present our conclusions in Section V. 



II. PRISONER'S DILEMMA WITH 
COOPERATION FACILITATORS: MODEL AND 
DYNAMICS 



In evolutionary game theory, two-player games can be 
interpreted as dilemmas of cooperation. In fact, the two 
possible strategies can be interpreted as "cooperation" 
(C) and "defection" (D). The paradigm of social dilemma 
is provided by the classic prisoner's dilemma (PD), whose 
main features are captured by the following payoff matrix 
giving the pairwise interaction between cooperators and 
defectors §, | |o) @: 



one has 



C 
D 



C D 

b — c — c 
b 



(1) 



where b and c respectively represent the benefit and the 
cost of cooperation, with b > c > 0. Here, without loss of 
generality, we assume that < b < 1. According to (Q), 
mutual cooperation leads to a payoff b — c > and mu- 
tual defection gives a payoff 0; whereas when one player 
defects and the other cooperates, the defector receives 
a payoff b and the cooperators gets — c. In the (classic) 
PD, the dilemma arises from the fact that each individual 
is better off not cooperating, even though mutual coop- 
eration enhances the population overall payoff. Hence, 
while cooperation is socially beneficial, defection is the 
only (strict) Nash equilibrium in the PD ^, || . 

In this work, we consider a finite population compris- 
ing N individuals on a complete graph (no spatial struc- 
ture). The number of cooperators and defectors is re- 
spectively denoted by j and k. In addition to coopera- 
tors and defectors, we consider that the population also 
comprises a fixed (small) number £ of "cooperation fa- 
cilitators" (I <C N). These facilitators cooperate with 
C— players and therefore enhance the reproductive poten- 
tial (fitness) of cooperators, while they leave the fitness of 
defectors unaltered, see below. Hence, while the number 
of cooperators and defectors in the population changes in 
time (j and k vary) , the total number of cooperators and 
defectors j + k = N — £ is conserved. According to the 
tenets of EGT, the variation in time of the number of co- 
operators and defectors depends on their average payoffs, 
7Tc and 7Tq respectively, obtained from the payoff matrix 
([j]). Here, since facilitators enhance ttc by cooperating 
with C individuals and have no (direct) influence on 7Tq, 



TTC = (b — C) — — ; C- 



iV-1 N-l 



7T D = b 



N-l' 



(2) 



where we have excluded self-interactions from the defi- 
nition of the payoffs 0. The population average pay- 
off is given by ff = (jVc + kito)/N. It is worth notic- 
ing that the expression of ttc now comprises a term 
(b — c)£/(N — 1) > reflecting the positive contribu- 
tion of facilitators to the cooperators payoff. In evolu- 
tionary dynamics, it is customary to add a baseline con- 
stant, here set to 1, to the payoffs tvc/d °f t nc spreading 
species |ll|, yielding the fitness of species C and D, 
respectively given by 



fc = 1 + 7T C = 1 + 6 



/ D = 1 + 7T D = 1 + b 



j + £- r(N - 1) - 1 
N-l 



N-l' 



(3) 



where, we have introduced the cost-to-benefit ratio r = 
c/b (with < r < 1) and have used k = N — j — t. Sim- 
ilarly, the average fitness of the entire population reads 
/ = 07c + kf D )/N = 1 + [6(1 - r)j - £]/N and grows 
linearly with the density x = j/N of cooperators. 

The size of the population being finite, the evolution- 
ary dynamics is modeled as a continuous-time birth- 
death process ^ |l7], 
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In this model, only pairs 
of cooperators and defectors interact (according to ([j])) 
and the stochastic dynamics is implemented as follows: 
(i) at each time step a pair of individuals is randomly 
chosen from the entire population; (ii) unless a pair of 
cooperator-defector is drawn, nothing happens; and (iii) 
if one picks a cooperator-defector pair, one of these in- 
dividuals is randomly chosen for reproduction (propor- 
tionally to its fitness) and the other is replaced by the 
newborn offspring. Hence, at each interaction the num- 
ber of cooperators increases or decreases by one. The 
time evolution of this birth-death process can therefore 
be described by the random variable j giving the number 
of cooperators and by the rates Tp associated with the 
transitions j — »• j ± 1, respectively. Here, we consider 



If = 



j(N-e-j) 

N(N - 1) 



^(/c/d), 



(4) 



where j(N — £ — j)/N(N — 1) accounts for the probabil- 
ity of picking a cooperator-defector pair, while are 
functions of the fitnesses (||) that encode the interactions 
(selection) according to the chosen "microscopic" update 
rule [pi. We here discuss the cases where <I' ± correspond 
to (i) the Fermi process (FP) fi~9| , po| and (ii) the Moran 
process (MP) |2| [12, 14, [2l| that are commonly used 
in EGT [§. 

Stochastic evolutionary dynamics and the influence of 
selection are generally characterized by the fixation prop- 
erties, namely the probability that a given species fixates 
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(takes over) the whole population and by the mean time 
for such an event to occur j|, |TJ[ [Tj|. In the absence of 
facilitators, fixation happens when only one species sur- 
vives and the population composition is uniform. Here, 
as the number of facilitators remains constant, fixation 
will be achieved when one of the absorbing states is 
reached and either all cooperators are replaced by de- 
fectors, or vice versa, resulting in a (non-uniform) popu- 
lation comprising £ facilitators and N ~ £ cooperators 
or defectors. In this work, we are particularly inter- 
ested in the probability <fij that, starting with j coop- 
erators, all defectors are eventually removed from the 
population and replaced by cooperators. As discussed 
in Sec. III. A., the fixation probability cj>^ is necessary to 
establish when selection favors cooperation replacing de- 
fection Ol . In the framework of the above birth-death 
process Wj, this probability obeys the backward master 
equation |2, 14, [L8) 



•27 



i + [1 - Tf - T+]tf 



with absorbing boundaries <pi 
formal solution of Eq. (|^) reads |l4| 



and <jy}f_ 



(5) 



1. The 



(2772? 



s-^N-t-l T-rr 
Z— m— 1 1 Lr, 



(0) 



Since the above birth-death process is a one-dimensional 
Markov chain, other quantities like the mean fixation 
times (MFTs) can, in principle, be obtained exactly, but 
yield unwieldy expressions When the popula- 

tion size N is large, it is often much more useful to de- 
scribe the fixation properties in terms of the diffusion 
approximation obtained in the continuum limit (TV 3> 1) 
by a second-order size-expansion of the master equation 
resulting in a Fokker-Planck equation [l5| [I?], [l8| . 
By denoting x = j/N and y = k/N the initial den- 
sity of cooperators and defectors, respectively; and with 
z = £/N being the fraction of facilitators in the pop- 
ulation, the (backward) Fokker-Planck equation (FPE) 
associated with (|J) reads |l7], 



Qback(x)(l) C (x) = 0, 



(7) 



where 



<!>{x) = <t>j /N 



and 



0back(aO = [T (x) - T (x)] 



dx 



(8) 



with T ± (x) = T.j N and, as usual, the density x changes 

by ±8 — ztN^ 1 at each cooperator-defector interaction. 
In the realm of the Fokker-Planck equation, the timescale 
is such that the time-step is S = N^ 1 . The formulation in 
terms of the FPE allows a neat connection with the mean 
field treatment of the dynamics. In fact, when N — > 00 
and all demographic fluctuations are negligible, the time 



variation of the density of cooperators is given by the 
drift term of (§1) IxtI , i.e. 



dx(t) 
dt 



T+{x) - T~(x) 



(9) 



x(l 



*)[* + (/c,/d)-*-(/c,/d)]. 



As for the classic PD, this rate equation admits two ab- 
sorbing fixed points, x — (no cooperators) and x = l—z 
(no defectors), but possesses no interior fixed point since 
# + (/c(z),/d(x)) + *-(fc{x),f D (x)) for the Fermi and 
Moran processes, see below. As discussed in what follows, 
the stability of these fixed points depends on the differ- 
ence between the cost-to-benefit ratio r and the fraction 
z of facilitators. 



III. DYNAMICS WITH THE FERMI PROCESS 

The stochastic dynamics of evolutionary games is of- 
ten conveniently modeled in terms of the so-called Fermi 
process (FP), see, for example, jl{|. In the FP, at each 
time-step two individuals are randomly drawn from the 
entire population and one of them reproduces at the ex- 
pense of the other that is replaced by the newborn off- 
spring. This happens with a probability proportional 
to the difference between the fitness of the interacting 
individuals and given by the Fermi function from sta- 
tistical physics. Since only the CD pairs interact, the 
dynamics with the FP is described by the birth-death 
process defined by_(@) and * ± = [1 + e T(/c-/o)]-i _ 



[l + eT^- 710 ']- 1 jTg|. With these expressions of * ± , one 
checks that * + (/c,/d) ^ *~(/c,/d) (since f c ^ /d, see 
Eqs. (|3|),([24[)), which confirms the absence of an interior 
fixed point in the mean field (continuum) limit. 



A. Fixation probability 

With (||) , the transition rates (f|) for the Fermi process 
read 



Tf 



1 



W-e-j) 

N(N-l) l + exp(±wjv)' 



with 

Vn = fo - fc = b 



(10) 



(11) 



pressure [|22| . 
< 1. In the 



This quantity measures the selection 
Clearly, — bz < vn < b(l — z) and \vn 
continuum limit N ^> 1, the densities x = j/N, z = £/N, 
and vn — > v = b(r — z) are treated as continuous quanti- 
ties, and the absence of self-interaction is ignored yielding 
the transition rates (|l0|) 



x(l 



x) 



1 + e 



(12) 
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FIG. 1: (Color online). Probabilities <fi and </>? for various 
z — i/N as function of j/(N — £) = x/(l — z), and evolu- 
tion with the Fermi process. Results of stochastic simulations 
(symbols) for (f>^ are compared with the predictions (curves) 
of (0) for z = (x, dot-dashed), 0.08 (□, solid), 0.12 (V, 
dashed). Similarly for (f>® with z = (o, thin dashed), 0.08 (o, 
solid gray), 0.12 (A, thin solid). The other parameters are 
N = 200,6 = 1.0, c = 0.1 (i.e. r = 0.1). Stochastic simula- 
tions are for the birth-death process defined by ( jl(i| ) and have 
been averaged over 2 x 10 samples. 



The rate equation corresponding to the mean field dy- 
namics with the Fermi process is obtained by using (|12j ) 
into (0), and is characterized by a single stable (absorb- 
ing) fixed point corresponding to a stationary density 



xq = 1 — z (no defectors) 
id = (no cooperators) 



if Vln ( 13 ) 
if v > v ' 



of cooperators. This means that cooperation prevails 
(x* = Xq, no defectors) in an infinitely large popula- 
tion comprising a fraction z of facilitators higher than 
the cost-to-benefit ratio r, i.e. when v < 0. However, as 
in the traditional PD, defection wins (x* — xq, no coop- 
erators) if z is less than r (v > 0). In other words, for 
cooperation to prevail (x* — xq) at mean field level, it 
is necessary that the density of facilitators compensates 
the cost of cooperation relative to its benefit. 

When the population size is finite, demographic fluctu- 
ations are important and the evolution is thus no longer 
aptly described by the mean field dynamics (^). In par- 
ticular, the mean field results (|l^) do not account for 
the nonzero probability that a single cooperator can in- 
vade and replace a (finite) population of defectors. Here, 
to investigate how the above mean field picture (p|),(|l3|) 
is altered by fluctuations arising in a finite population, 
we compute the probability that defection is eventu- 
ally replaced by cooperation in a population comprising 
initially j cooperators and N — £ — j defectors. Since 



T+/T7 



this probability can be obtained explic- 



itly using (^J) and, when vn =/= j22j, one finds 



e N(l-z)v N _ i ' 



(14) 



while fixation probability of defectors is simply given by 

= 1 - <jyj = ( e ^(i-*>N - e Jv N y( e N(i-z)v N _ l y As 

shown in Fig. [j], where results of stochastic simulations 
obtained using the Gillespie algorithm [^3| are reported, 
the predictions of (14) are in excellent agreement with 
numerical simulations. The expression (flj) implies that 



1 -e 

(e> v - 



-3\vn\ 

1) e" 



if v N < 

N(l~z)v N if WAr >0, 



(15) 



when -/V|«at| ;g> 1. In particular, the cooperation fixation 
probability starting with a single cooperator reads 



1 



z N(l-z)v N _ i 



1 



-\vf 



(g"iv _ \\ e -N(l~z)vp, 



i ! UjV< n (16) 
if v N > 0. v ; 



The findings (|l4l)-(|l6|), summarized in Fig. [j], illustrate 
how a small fraction z of cooperation facilitators affects 
the fixation probabilities in a large, yet finite, popula- 
tion with an initial density of cooperators comparable to, 
or larger than, the density of defectors: When vn < 
(z > r), the fixation probability of cooperators is much 
higher than that of defectors, <f>$ 3> <ffj , and the spread 
of cooperation is thus efficiently promoted by facilitators. 
Yet, it is worth noticing that defectors still have finite 
probability to fixate even when z > r (and x <C 1 — z), 
contrary to the mean field prediction (13). The opposite 
situation arises when vjy > (z < r), as shown in Fig. [I]. 

The results ([l5]) and (|l6|) can also be used to assess the 
influence of selection on the evolutionary dynamics [^): 
Following the seminal work of Ref. jF3|, we can estab- 
lish when selection favors cooperation (C) invading and 
replacing defection (D). Selection is said to favor the re- 
placement of D by C if the fixation probability <f>i of a 
single cooperator in a population of N — i — 1 defectors 
is greater than in absence of selection pressure (vn = 0) 
when (ft tVri=0 = (N - I)- 1 |§. With <JX6|) , this yields 
the condition 1 — e - '"™' > (N(l — z))^ 1 that is generally 
satisfied in large populations under non-vanishing selec- 
tion pressure. An interesting result arises when the se- 
lection intensity is weak and the population size is large, 
i.e. \vn\ -> M <C 1 and N 3> 1. In such a limit, 
cf>\ ~ \v\ when z > r (see Fig. || where \v\ = 0.02) 
and selection favors cooperation replacing defection pro- 
vided that b(z — r) > (N(l — z))^ 1 . Moreover, selec- 
tion favors C invading D when fc > /d [jl3). With (||), 
this yields the condition z — r > (1 — r)/N. There- 
fore, under weak selection (|u| <C 1 and N » 1) se- 
lection favors the invasion and replacement of D by C if 

Since < b < 1, one 



r > 



max ( 1 



6(1- 



N ' 

has 1 — r < (6(1 — z))^ 1 and cooperation invading and 
replacing defection is favored by selection provided that 



b(z-r)(l-z) > N~ 



(17) 
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One can also use the results (|l6|) to determine the cir- 
cumstances under which defection is evolutionary stable. 
In fact according to p3[ , and as natural extension of 
the concept of evolutionary stability for infinitely large 
populations and deterministic evolutionary dynamics [0 , 
D is evolutionary stable in a finite population if (i) se- 
lection opposes C invading D, implying fc < fo, he. 
z — r < (1 — r)/N; and if (ii) selection opposes C replac- 
ing D, i.e. (f>i < (A — I) -1 . The condition (ii) is clearly 
always satisfied when \vjy\ is finite, and in this case de- 
fection is evolutionary stable if z — r < (1 — r)/N. In the 
weak selection limit where \v\ <C 1 (with N ^> 1), the 
condition (ii) yields z — r< (bN(l — -z)) -1 . Hence, defec- 
tion is evolutionary stable under weak selection in a large 

population if z — r < min ^1 — r, b ^i_ z ^ ^ = (1 — r)/N . 
Since r < 1, this clearly means that defection is evolu- 
tionary stable and is the dominating strategy when z < r. 
It is worth noticing that in the limit of an infinite popu- 
lation, A — > oo, one recovers the mean field results (|13|): 
cooperation prevails only if z > r, according to ([T^), and 
defection dominates otherwise. 

The meaning of the results (|l5|)-(|l7|) is illustrated 
in Fig. U where has been computed in populations 
comprising a small initial number of cooperators (j — 
1, 10) and an excellent agreement with ( |l5| ) and ( |l6|) 
has been found. In Fig. || \vn\ <C 1 and we notice that <^ 
increases linearly in x — j/N <C 1, with a slope steeper 
than (1 — z)~ x when z > r and selection favors coopera- 
tion replacing defection. The slope is less than (1 — z) _1 
when z < r and the fixation of cooperation is opposed 
by selection. To further appreciate the implications of 
(|l4|)-(|l7|), it is useful to compare these results with those 
obtained in the absence of facilitators. Putting z = 
in (|l4|)-(|l7|), one recovers the results for the classic PD 
when cooperation fixation probability vanishes exponen- 
tially with the population size A: 4>~ z=0 ~ e -( N -J) c anc [ 
4>l z= v~e- Nc |,|H (see Fig. 0). " 

Our findings therefore demonstrate that facilitators 
greatly influence the probability that cooperation pre- 
vails and are summarized in Fig. [l]. As illustrated in 
that figure, the influence of facilitators crucially depends 
on the difference between their density z and the cost-to- 
benefit ratio r: 

i. When vm < 0, the fixation of cooperators is likely 
(but not certain) even when they are initially in 
minority, i.e. even when initially x = j/N < 1/2. 

ii. When vn < and A|uiv| S> 1, the fixation prob- 
ability of a single cooperator is generally higher 
than in the absence of selection pressure (vn = 0). 
In this case, with z > r, selection favors coopera- 
tion invading and replacing defection, see ( |l6| ) and 
Figs. [l] and |2[ Furthermore, under weak selection 
pressure and in a large population (\v\ <C 1 and 
A 1), the fixation probability of a single cooper- 
ator is independent of A, <f>^ ~ z — r. In this case 
invasion and replacement of defection by coopera- 
tion is favored by selection if (p7|) is satisfied. 




j/N 



FIG. 2: (Color online). Probability <j£ as function of j/N 
when the initial number of cooperators is j = 1 — 10 with 
N = 500, and j = 1 - 5 with N = 200. The dynamics is im- 
plemented according to the Fermi Process with (p"(i|). The re- 
sults of stochastic simulations (symbols, averaged over 2 x 10 
samples) are compared with (114) (curves/lines). Parameters 
are b = 1.0, c = 0.1 (i.e. r = OA), and (N, z) = (500, 0.12) (o), 
(200,0.08) (o). Here, ^ ~ 0.0182 (o) and $ ~ 2.88 x 
10 -4 (o), see text. The dashed/dotted/thin lines correspond 
to (j>f ~ jv N /(e N(1 ~ z)VN - 1) with (N,z) = (500,0.12) 
(dashed) and <f>f, VN =o = j/{N-€) for (N,£) = (200, 16) (thin) 
and (N,£) = (500,60) (dashed-dotted). 

iii. When vn > 0, selection always opposes cooper- 
ation replacing defection. In this case, while de- 
fection is evolutionary stable and is the dominat- 
ing strategy when z < r, the cooperation fixation 
probability is exponentially enhanced by a small 
fraction of facilitators. Yet, cooperation is likely to 
fixate only if defectors are initially outnumbered by 
cooperators, i.e. if j ^> k, as illustrated in Fig. [I]. 



B. Mean fixation times 

Another quantity of great interest to unveil the influ- 
ence of facilitators in the evolutionary dynamics of the 
PD is the (unconditional) mean fixation time. This quan- 
tity gives the average time necessary to reach one of the 
absorbing boundaries, i.e. a population composition with 
either or N — i cooperators. The unconditional mean 
fixation time (MFT) , Tj , for a system comprising initially 
j cooperators obeys the following backward master equa- 
tion g |lj, [l8) (where the time-step is S = A -1 ) 

Ti=* + Trr^! + T+ Tj+1 + [1 - Tr - T+]t 3 , (18) 

with boundary conditions tq — r^v-f = 0. In principle, 
this equation can be solved exactly but the final result is 
cumbersome and not enlightening. Here, in the contin- 
uum limit JV> 1, we work with the continuous quanti- 
ties x — j/N, z — l/N and v — b(r — z), and adopt the 
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approach of the diffusion theory [12 uM. The diffusion 



approximation is known to be particularly suited to ana- 
lyze the dynamics under weak selection, which here cor- 
responds to the regime where |v| <C 1 |, ||. Exact 
methods (when available) or other approximations [ [l4| , 
e.g. the WKB approach PQl ) are particularly useful to 
deal with the case of strong selection intensity and/or 
with phenomena like metastability. In the realm of the 
diffusion theory, the transition rates of the FP are given 
by and the fixation probability of cooperation is ob- 
tained by solving (0) which yields 



ctf{x) 



1 



e N(l-z)v _ 1 ' 



while for defection the probability is (jP{x) = 1 — <j> (x). 

Similarly, the unconditional MFT is obtained by solv- 
ing the backward FPE Gback(x)T(x) = -1 [00, i.e. 



x(l 



tanh 



d 
dx 



J_d*_ 

2N dx 2 



t(x) = 1, (19) 



with the absorbing boundary conditions r(0) = r(l — 
z) = 0. When the drift and diffusive terms are of the 
same order, i.e. when \v\ ~ AT -1 <C 1, it follows from 
Eq. (ph that the MFT scales linearly with N, i.e. 



-(x) = NF v {x). 



(20) 



The scaling function can be obtained explicitly by solving 
Eq. (19) using standard methods, see e.g. pTj . For in- 
stance, when the initial density of cooperators and defec- 
tors is the same, x = y = (1 — z)/2, and \v\ ~ N^ 1 <C 1, 
one finds 



1 



-(2+z)g 



3 (2+z)9 



q(l - z){l + e-i^- z >, 
x { 7E -ln2-Ei(-(l-z) g )} + e( 1+2z ) 9 
x {Ei((l 7e - In (2g)} (21) 
+ e 3zq {Ei((l - z)q) - Ei(2(l - 
+ e 3 " {Ei(-2(l-z) (? )-Ei(-(l-z) (? )} 

+ e (1+z)<? (e 9 -e z<? ) ln(l-z) 



J — OO U 

= 0.5772... 



where q = N\ tanh (u/2)| ~ AT|u|/2, Ei(x) ee 
denotes the usual exponential integral and 7e 
is Euler-Mascheroni constant. While the expression ofT v 
is usually cumbersome, some useful properties can be di- 
rectly inferred from ([H]). In fact, as Eq. ( |T9| ) is invariant 
under the transformation (x,r) — > (1 — z — x, 2z — r), 
one has !F v {x) — J-L^l — z — x) when z is kept fixed. 
The unconditional MFT in the Fermi process is there- 
fore characterized by the symmetry 



"(x)|r =t(1 — z — x)\, 



(22) 



where, on the right-hand-side r is replaced by r' = 2z — r 
and v transformed into — v, with z kept fixed. Further- 
more, when r = c/b is kept fixed but z varies, ( |l9| ) is 




0.4 0.6 

j/(N-£) 

FIG. 3: (Color online). Mean fixation times as function of 
j/(N — £) — x/(l—z) for the evolution with the Fermi process. 
Results of stochastic simulations (symbols) for r are compared 
with the solution (curves) of Eq. (M) for z = (o), 0.08 (o, 
solid black), 0.12 (x, solid gray). We also report the numerical 
results for the conditional MFTs r c for z = 0.12 (□) and r D 
with z = l/N = (V), 0.08 (A). The other parameters are 
N — 500, b = 1.0, c = 0.1, i.e. r = 0.1. Stochastic simulations 
are for the FP with rates ( |l(i| ) and have been averaged over 
2 x 10 5 samples. 



invariant under the transformation 



= 2r — z and 



x — > 1 — z' — x, while the boundary conditions become 
t(1 - z') = and t(z - z') = r(-2v/b) = 0. In the weak 
selection regime \v\/b = \z — r\ <C 1, the second boundary 
condition can be approximated by t(z — z') ~ r(0) = 0, 
which allows a mapping onto (|l9|) that yields: 



t{x)\ z ~ r(l - z' - x)\ 2 



(23) 



with r = c/ b fixed. The comparison between the so- 
lution of (119) and the results of stochastic simulations 
(for the FP with rates ©) reported in Fig. | shows 
that the diffusion approximation aptly captures the func- 
tional dependence of r, even though some deviations (of 
about 10%) can be noticed. These deviations stem from 
the self-interaction terms that are excluded from jl0| ) 
but not in the continuum limit ( |l2| ) [e.g. in Fig. ^ one 
has v N ^ -0.0182 and v = -0.02 when z = 0.12, and 
v N ~ 0.0218 and v = 0.02 for z = 0.08]. More im- 
portantly, the scaling ( pp| ) and the relationship (^) are 
confirmed by the numerical simulations of Fig. ^. In fact, 
in Fig. U we notice that t(x) is a humped function with a 
maximum well separated from the absorbing boundaries 
and located at x/(l — z) < 1/2 when z > r and, while 
r scales linearly with N, the presence of facilitators in- 
creases the unconditional MFT and its maximum value 
at the hump. 

In addition to the unconditional MFT, it is also rele- 
vant to consider the mean time to specifically reach one of 
the absorbing boundaries. Hence, the conditional mean 
fixation times t c (x) and t d (x) respectively give the av- 
erage time to reach the absorbing boundaries x = 1 — z 
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and x = @ f(|. As for the unconditional MFT, 
these quantities can be obtained from a backward FPE in 
the realm of the diffusion approximation. In fact, t c (x) 
obeys Gback(x)[(l) C (x)t c (x)] = —(j) C (x), with the absorb- 
ing boundaries ^ c (l - z)r c (l - z) = c (O)r c (O) = @. 
Since <j> D {x) = 1 - <p c {x) and, from @, <p D {x) = 
4> c (l — (2r — z) — a;), the conditional MFTs in the regime 
\v\/b -c 1 (weak selection) are related by the relation- 
ship t c (1 + z - 2r - x)| z ~ T D (x)J 2 _ >2 ' = 2r-2 where r 
is kept fixed, as illustrated in Fig. |. Furthermore, one 
has </> D (cc) ~ 1 when x — > and u > (z < r), while 
</> c (:r) ~ 1 when x — > 1 and u < (z > r). This im- 



-C/ 



r"(x) when v > and a; — >• 



As 



plies that t(x) '— ~. , , , 

I r (x) when w < and x — > 1 

shown in Fig. ||, r c (a;) decreases while t d {x) increases 
monotonically with x/{\ — z). 

The influence of facilitators on the unconditional and 
conditional MFTs is summarized in Fig. ||. We have 
found that in the PD with cooperation facilitators all 
conditional and unconditional MFTs scale linearly with 
the population size N when \v\ ~ iV" 1 (weak selection). 
While a similar scaling is also obtained in the absence of 
facilitators, the MFTs at a fixed value x/(l — z) are found 
to be significantly increased by the presence of facilita- 
tors. Hence, the presence of cooperation facilitators has 
the quantitative effect to prolong the coexistence and the 
competition between cooperators and defectors before an 
absorbing state is reached, see Fig. [|. 



IV. DYNAMICS WITH THE 
FITNESS-DEPENDENT MORAN PROCESS 



The stochastic dynamics of evolutionary games is of- 
ten implemented in terms of the Moran process, see 
e.g. 1^1 , that was originally introduced in population 



genetics||12j, |2l| . In its essence, the Moran model is a 
birth-death process where one randomly picked individ- 
ual produces an offspring proportionally to its fitness rel- 
ative to the population average fitness. The resulting 
offspring then replaces another individual that is ran- 
domly picked to be removed from the population whose 
size is therefore conserved. Here, as the interactions are 
between cooperators and defectors, the Moran process is 
implemented with ^ + — fc/f an d = .fn/.f in (^). 
Since /c ^= fo when v ^ [see (24) and [P2| 1, one veri- 
fies that v I' + (/c,/d) 7^ ^"(/c/d) implying the absence 
of an interior fixed point in the mean field (continuum) 
limit. The Moran process is usually investigated when 
the selection intensity is weak, both for technical con- 
venience (the mathematical treatment simplifies gre atly) 
and for the biological relevance of such a limit |12|, [l3| . 
In this section, the stochastic dynamics with the Moran 
process is investigated in the weak selection limit, where 
|u| = 6|r — z| <C 1, using the diffusion approximation. 



A. Fixation probability 

In the continuum limit, the fitnesses (|^) become 

f c (x) = l-v + bx and f D (x) = 1 + bx, (24) 

with f(x) — I — z + 6(1 — r)x. The transition rates for 
the Moran process thus read 



T+/- 



(x) =x{l-z-x) - . 

/0) 



(25) 



With ( |24| ) and (|25|) , the mean field dynamics is described 
by the rate equation (0) whose properties are similar to 
those discussed for the Fermi process. In particular, the 
rate equation (Q) for the Moran process is also character- 
ized by a single stable (absorbing) fixed point x* = xq 
(no defectors) if v < and x* — xq (no cooperators) if 
v > [see @]. 

To understand how the combined effect of nonlinear 
selection and demographic noise alters the mean field de- 
scription, we now compute the cooperation fixation prob- 
ability in the realm of the diffusion approximation. In 
such a setting, the fixation probability <fi c (x) is given by 
the FPE (0|f) with the boundary conditions </> c (0) = 
and 4> c (l — z) = 1. The solution of (Q) is given by [|l7| 



f{x) 



where, with (g4j) and (|5| 



Jo du x(u) 
Jo~ Z du x(u) 



(26) 



X(«) = exp J dsS. 



exp 2Nv 



T+(s)-T-(s) 
T+(s)+T-(s) 
ds 



2bs + 2 - v 



(27) 



Introducing (|27]) into ([26|) and performing the integrals, 
one obtains 



4> c (x) 



l+Nv/b 



1 



l+Nv/b 



(28) 



As shown in Fig. this result is in excellent agreement 
with numerical simulations and exhibits the same quali- 
tative features obtained for the Fermi process (compare 
with Fig. |l|). The finding ( p8| ) implies that in the weak 
selection limit where \v\ <C 1 and N\v\ S> 1, one has 



<t> c {x) 



1 - (I + bx)- N ^- r ^ ifz> 



V l+fc(l-i+I 



N(r-z) 



if r > z. 



(29) 



In particular, the probability that cooperation fixates 
starting with a single cooperator, when z > r is given by 
\iuiN X ^i 4> c (x) = 1 — e ~ \v\. We therefore recover the 
result derived from (JL6h for the Fermi process. Clearly, 
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1.0 



FIG. 4: (Color online). Probabilities (j) C {x) and <jP( x ) f° r 
various 2 as functions of x/(l — z), and the dynamics with 
the Moran process (p5|). Results of stochastic simulations 
(averaged over 2 x 10° samples) are compared with (0) for 
z = (x, dot-dashed), 0.08 (□, solid), 0.12 (V, dashed). 
Similarly for <jP(x) with a = (o, thin dashed), 0.08 (o, solid 
gray), 0.12 (A, thin solid). The other parameters are N = 
200, b = 1.0, c= 0.1 (H = 0.02). 



FIG. 5: (Color online). Rescaled mean fixation times as 
function of x/(l — z) for the evolution with the Moran pro- 
cess (|l^). Results of stochastic simulations for r are com- 
pared with ( |3C| ) for z = (o, dashed), 0.16 (o, solid black 
curve), 0.24 (x, solid gray). Numerical results for the con- 
ditional MFTs r c with z = (*), 0.16 (D),0.24 (+) and for 
r D with z = (V), 0.16 (A), 0.24 (•). The parameters are 
N = 200, b = 0.25, c = 0.05 (r = 0.2 and |«| = 0.04). Stochas- 
tic simulations have been averaged over 2 x 10 samples 



this implies that under weak selection the fixation of a 
single cooperator is favored by selection if the non-trivial 
condition ([171) is satisfied. Again, it is instructive to com- 
pare (E8h , (|29|) with the result obtained in the absence of 



facilitators, when </> (x)| z= o 



l+b(x- 
1+6(1- 



r/2) 
7751 



A'r 



decays to 



zero exponentially with AT. The influence of the facilita- 
tors on the fixation probabilities for the Moran process 
is summarized in Fig. ||, where the same features as in 
Fig. ^ are recognized and summarized as follows: 

i. The fixation of cooperators is likely (but not cer- 
tain) when the density of facilitators is higher than 
the cost-to-bencfit ratio (z > r). 

ii. When \v\ <C 1 and N\v\ 3> 1, selection favors coop- 
eration invading and replacing defection if (17) is 
satisfied. In particular, the fixation probability of 
a single cooperator is liniAra;- 



iii. When z < r, selection opposes cooperation replac- 
ing defection but the fixation probability of coop- 
erators is exponentially enhanced by the presence 
of facilitators. 



B. Mean fixation times 

In the realm of the diffusion approximation, the uncon- 
ditional mean fixation time r obeys the backward FPE 
5back(x)r(x) = —1, with the absorbing boundary condi- 
tions t(0) = t(1 — z) = 0. In the weak selection regime 



c<()Cl and continuum limit, with (25), one has 



T + {x)-T-{x) 
T+(x)+T-(x) 



■ x(l — z — x), 



1 



z 



x). 



With these expression, the backward FPE for the uncon- 
ditional MFT reads 



x(l — z — x) 



dx 



JVtf 



t{x) = -1, (30) 



with r(0) = t(1 - z) = 0. When \v\,b < 1, Eq. @ 
coincides with the FPE ( |i~9| ) for the Fermi process with 
an effective population size N(l — z)/2. The solution 
to (30) can therefore readily be obtained from (19) and 
@. In particular, we infer from (p3) that the MFT 
scales linearly with N(\ — z)/2 when \v\ ~ N^ 1 , yielding 



t{x) 



N(l - z) 



Fvix), 



(31) 



where J- V {x) is the scaling function ( p0[ ) obtained for the 
Fermi process. This function still satisfies the symme- 
try J- V (x) = J-- V (l — z — x) yielding r(x)\ r = r(l — z — 
x)| T ._ >r / = 2 2 _ r , when z is kept fixed, as in the Fermi pro- 
cess. In the same manner, from (^l]) and (p3|), we infer 



t(x) 



1-z 



1 + z - 2r 



r(l 



=2r-z 



(32) 



when r — c/b is kept fixed and z is transformed into 
z 1 = 2r — z. The solution of (|30|), as well as the rela- 
tionships (|3l|) and (|32|), are in excellent agreement with 
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the results of stochastic simulations reported in Fig. |5j. 
As for the FP, we can also consider the conditional mean 
fixation times and it follows from (|3l]) and ( p3| ) that for 
the Moran process the conditional MFTs are related by 

( ld Tfr I ) r c (l + z ~ 2r ~ X )U - T D {x)\ z ^ z , =2r - z where 
r is kept fixed, in agreement with the results of Fig. [|. 

The influence of facilitators on the MFTs with the 
Moran process is summarized in Fig. ||, where the MFTs 
rescaled by a factor (N(l — z)/2)~ 1 reproduce the same 
qualitative behavior obtained for the Fermi process (com- 
pare with Fig. |3|) and t(x) is a humped function with a 
pronounced maximum. Again, all MFTs scale linearly 
with N (in the weak selection limit). Yet, the compar- 
ison with the results for z = reveals that, at a fixed 
value of x/(l — z), the presence of facilitators increases 
the MFTs, see Fig. hi Also, we notice that the monotonic 
dependence of r c and r D on x is essentially independent 
of the sign of v ^ (in Fig. [|, v = ±0.04 and v = 0). 

V. SUMMARY AND CONCLUSION 

In this work, we have proposed and investigated an al- 
ternative scenario leading to the spread of cooperation 
in social dilemmas. We have considered the evolution- 
ary dynamics of the prisoner's dilemma (PD) game in 
the presence of a small number of cooperation facilita- 
tors. These individuals participate in the dynamics only 
by enhancing the fitness of cooperators. The influence of 
facilitators on the evolutionary dynamics has been char- 
acterized by computing the model's fixation properties 
in a finite population of size N. Here, fixation occurs 
either in the state with only defectors (as in the classic 
PD), or in the state where the entire population is com- 
prised of cooperators and facilitators. The dynamics has 
been implemented with the Fermi and Moran processes 
and the same qualitative results have been found, which 
demonstrates the robustness of our findings. Our ana- 
lytical approach, corroborated by stochastic simulations, 



is based on an exact treatment and on the diffusion ap- 
proximation (Fokker-Planck equation) of the underlying 
birth-death process. 

Our main results concern the fixation probabilities, 
whose properties crucially depend on whether the frac- 
tion of facilitators z is more or less than the game's cost- 
to-benefit ratio r. When z > r, we have shown that 
facilitators are very efficient in promoting the spread of 
cooperators whose fixation is likely (but not certain, con- 
trary to the mean field predictions) in a large population 
with comparable initial densities of defectors and cooper- 
ators. Furthermore, when the selection intensity is weak 
and JV > 1, we have demonstrated that the invasion and 
replacement of defectors by a single cooperator is favored 
by facilitators and selection if b(z — r)(l — z) > N^ 1 
(where < b < I is the cooperation payoff benefit). 
When z < r, defection is evolutionary stable and is the 
dominating strategy. In this case, while cooperation is 
unlikely to fixate, the fixation probability of cooperators 
is still exponentially enhanced by the presence of facili- 
tators. We have also studied the (unconditional and con- 
ditional) mean fixation times in the weak selection limit 
and found that these quantities grow linearly with the 
population size. While a similar scaling is also obtained 
in the absence of facilitators, their presence has the ef- 
fect of significantly increasing all the mean fixation times 
and hence to prolong the coexistence of cooperators and 
defectors. 

In conclusion, this work demonstrates that the pres- 
ence of a small number of cooperation facilitators can 
effectively enhance the spread of cooperation in a simple 
model of social dilemmas and prolong the coexistence of 
competing species. The influence of facilitators is partic- 
ularly drastic when their abundance exceeds the game's 
cost-to-benefit ratio, in which case cooperation is gener- 
ally the strategy favored by selection in large populations. 
These findings pave the way to further investigations of 
the influence of facilitators in other social dilemmas, e.g. 
with mixed strategies and/or in spatial settings. 
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